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Abstract. Based on liep-ph/0510121 , we discuss further the numerical study 
of classical SU(2) 3+1-D Yang-Mills equations for matter produced in a high 
energy heavy ion collision. The growth of the amplitude of fluctuations as 
exp (F-y/ g^i^r) (where jjL is a scale arising from the saturation of gluons in 
the nuclear wavefunction) is shown to be robust over a wide range of initial 
amplitudes that violate boost invariance. We argue that this growth is due 
to a non-Abelian Weibel instability, the scale of which is set by a dynamically 
generated plasmon mass. We discuss the relation of P to the prediction from 
kinetic theory. 
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1. Introduction 

One objective of the experiments done at ultra-relativistic heavy-ion colliders such 
as RHIC and, in future, the LHC, is to understand the properties of very hot 
and dense partonic matter in QCD. This requires understanding how the coherent 
wavefunctions of the incoming nuclei decohere, possibly forming a thermal Quark 
Gluon Plasma (QGP). At high energies, the small x (or wee) partons determine 
the properties of nuclear wavefunctions. Their properties can be formulated in 
an effective field theory called the Color Glass Condensate (CGC) A semi- 

hard scale Qs{x) » Aqcd, the "saturation" scale [EI, arises naturally in this 
approach and grows with energy; weak coupling techniques are therefore feasible. 
Furthermore , the small x wavefunctions of the incoming nuclei can be treated as 
classical fields with large occupation numbers [ (Sj . This enables the description of 
nuclear collisions in terms of solutions of classical Yang-Mills equations with the 
fields representing the small x partons, and light cone currents describing the hard 
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valence partons [ 4 . The latter can be modeled as 

= S''+pi{x^)S{x-) + 6''-p2{x^)d{x+), (1) 

where the color charge densities pi_2 of the two nuclei arc independent sources 
of color charge {x^ = {t ± z)/2). The 5- function sources ensure that the fields 
are boost invariant, namely, independent of the space time rapidity 77, defined as 
77 — atanh(z/t). The Yang-Mills equations can therefore be expressed in terms of 
the two transverse directions (xj^) and the proper time r, defined as r = \/t^ — . 
The initial conditions can be determined by matching the Yang-Mills equation in 
the four light cone regions, at t = 0, to determine self-consistently the fields in the 
forward light cone in terms of those before the collision. These latter classical fields 
can be computed analytically in the CGC framework. 

In the McLerran- Venugopalan model (MV) [El for large nuclei, the sources of 
color charge are Gaussian distributed 

{p1{^^_)p]{Y^)) = g^pH,,5'^'5\^^ - yx) , (2) 

where p is the dimensionful momentum scale in the problem. This scale is closely 
related to the saturation scale Qs which, in the classical effective theory, is defined 
as Ql = g^p?Nch\{g'^p/h.QCT))l'2-Tr- For the initial conditions corresponding to 
configurations of color sources of each of the two nuclei, the Yang-Mills equations 
can be solved numerically, and the final gauge field configurations averaged over the 
sources, to determine energy and number distributions [Ej- 

It is clear that for all finite \/s the ansatz of 5-function sources in Eq. has 
to be modified in order to implement the fact that the nuclei will not be contracted 
into infinitely thin sheets. More important, however, are the effects of quantum 
corrections, which may be of order unity over rapidity scales ~ I/as- These are not 
included in the MV model but arise from small x quantum evolution of the classical 
fields Consequently, one must deal with functions having a finite width in a;^, 
respectively. For a single nucleus it is still possible to solve the Yang-Mills equations 
classically and obtain the Weizsacker- Williams fields. For two nuclei, however, the 
problem becomes more involved, simply because the nuclei will interact for a finite 
time and the single nucleus solutions before the collision will be distorted during 
this time span. Ignoring the details of this process, the main difference with respect 
to the cases considered so far [ 5 will be the emergence of rapidity fluctuations and 
consequently a breaking of the boost-invariance of the small x fields. In what follows, 
we will concentrate on studying the effect of rapidity fluctuations by numerically 
solving the Yang-Mills equations after the collision (based on Ref.[|Sj); to keep 
the analysis as simple as possible, we assume the initial distortions of exact boost- 
invariance to be very small. The reasons for this are two-fold. One is to connect 
our results to published results [ 5 . The other reason, as we will discuss, is to 
study the effects of the Weibel instabilities over several decades in the magnitudes 
of amplitudes. 

This work is organized as follows: In section |21 we discuss why this study is 
connected to the phenomenon of non- Abehan plasma instabilities ([□ El UHl CH 
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ll^p. before providing details of our setup in section 13 Our results are presented in 
section 0] 

2. Motivation 

At earliest times rQg < 1, typical gluon occupation numbers are large and thus the 
system is described most appropriately in terms of nonlinear gluonic fields, which 
should be accessible by simulating classical Yang-Mills dynamics [Ej- However, 
because of the rapid longitudinal expansion, the gluon occupation number drops 
until the non-linearities become so weak that the hard modes (p ^ Qs) can be 
described as on-shell particles. In this regime (which should be reached for tQs > 1), 
the dynamics of the system is in terms of hard particles coupled to soft (fc <C gQs) 
fields, so a Vlasov-type kinetic approach should be appropriate to describe the 
system. Consequently, at times tQs 1, one would expect both classical and 
kinetic theory descriptions to offer a fair approximation of the system dynamics. 

Another consequence of the longitudinal expansion is that the typical longitu- 
dinal gluon momentum, for a fixed slice in rapidity, becomes smaller as 1/t. 
Since the transverse gluon momentum stays approximately constant p± ^ Qs, the 
gluon distribution function /(p) tends to become more and more anisotropic (until 
scatterings become important at very late times). 

Using a Vlasov approach it has been shown [ 13 El IIHI 1111 IT^ that when 
expansion effects are negligible, systems with an anisotropic momentum-space dis- 
tribution function are subject to the presence of so-called plasma-instabilities, with 
a typical exponential growth rate 7 proportional to the soft scale moo. 



in the limit of very strong anisotropics [ 9 . These instabilities manifest themselves 
as exponentially growing magnetic fields which in turn reduce the momentum- 
anisotropy, both by transferring energy from hard to soft excitations as well as 
by bending hard particle trajectories. 

Because of the relation between classical field dynamics and kinetic theory at 
times tQs ''^ 1 conjectured above, one would expect to see some manifestation of 
these instabilities when simulating classical Yang-Mills dynamics. To simulate an 
expanding metric, we solve the Yang-Mills equations in (r, 77,x_l) co-ordinates. In 
momentum space, the conjugate momenta are (fc"^, fc^, kj^), respectively. The Yang- 
Mills fields for "soft" fc,, modes will thus be sensitive to anisotropic distributions 
of modes in (kj^,/c^), thereby triggering an instability of the Weibel type. It was 
predicted by Arnold, Lenaghan and Moore [ 9 that in an expanding system such 
an instability would grow as exp (y^) rather than exp(T). As shown in Rcf. 
this is precisely what happens. Below, we discuss in some detail the setup of the 
numerical problem and some of the results. 
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3. Setup 



In yl"^ = gauge, the gluonic part of the QCD action has the form 



S 



dTdr]dx_\_TTT 



F2 

TTj 



F2. 



= I dTdr]dx±C, 



(4) 



where F^^i, = d^^Ai, — d^A^^ + ig[Afj^, A^,] is the field strength in the fundamental 
representation with F^'' = F^^Ta and [Ta,n] = ifabcTc, Tr t'^t'' = In the 

following we shall ignore effects of the current j^. This is justified if we limit 
ourselves to a small region around 77 = 0. With this restriction in mind we derive 
the conjugate momenta from the Lagrangian Eq.l^J, 



E, 



dC 



TdrA,, Er, 



dC 



d{drA,) " " d{drAn) 

with which we construct the Hamiltonian density 

'e} 



n = E,{drA,) + E^idrAr,) - C = Tv 



F2. 

r]i 



te: + tf: 



(5) 



(6) 



Here transverse coordinates x, y have been collectively labeled by the Latin index i. 
Using finally dH/dE^ = drAf^, dTL/dA^ = —drE^ , Hamilton's equations for the 
fields and their conjugate momenta are 



drA., 



E., 



drA,^ = rErf, drEi = tDjF, 



T-^D„F 



Since it will be used in the following, we also introduce two relevant components 
of the stress-energy tensor in (r, x, y, rf) coordinates. 



_^ j,yy ^ 2r Tr [Fly + E^] 
^2j.,v ^ r-'T:r[F^, + Ef]-rTv[F, 



2 

xy 



Ei 



(7) 
(8) 



3.1. Initial conditions - Boost-Invariant Case 



In the case of exact boost invariance, one obtains the initial conditions by matching 
the equations of motions before the collision (when there are only undisturbed single 
nucleus solutions) at the point a;* = and along the boundaries a;"*" = 0, a;~ > 
and x~ = 0, x"*" > 0. Omitting the details worked out in the result for the 

fields and momenta at time r = is 

Ai{xi_) = Q;i,i(a:;j_) + a2Aix±) A,j{x±_) = 
£iix±) = £ri{xj_) = ig[ai,i{x±), a^xj.)], 

"1,2 = -gUiad'Ul^, U,,2 = 7'exp(-z/^^rfa:'±Ai,2(xi,x'±)) (9) 

where V denoting path-ordering, AAi_2(a;_L, x^) = —pi,2{x±)S{x^) and pi^2 are to 
be determined from Eq.lj^J. 
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3.2. Initial conditions including rapidity fluctuations 

Ignoring the details of the initial rapidity profile we simply start with the boost- 
invariant field configuration and disturb it by adding small random rapidity vari- 
ations. Specifically, one has = Ai, A^ = 0, = SEi, E^j — £^ + 5Ej^, with 
DiSEi + DfjEjj = at initial time t = Tinit- The rapidity dependent functions SEi, 
SEjj are constructed as follows: for SEi we draw random configurations SEi{x±) in 
the transverse plane, < 6Ei{x±)SEi(y±) >= S'^{x± ~yi_) and subsequently multiply 
by a random function f{'q) = drfF^rj) with dimensionless amplitude A <C 1, 



to get SEi{xi, rf) = f{r])6Ei{xi); 5E^ is then constructed as 5Er^ = —F{f])DiSEi{xi). 

Thus, by construction, one has added random rapidity fluctuations of amplitude 
A to the system which obey Gauss's law. An advantage of this construction is of 
course that one can use periodic boundary conditions in the 77 direction, which for 
a lattice simulation is somewhat more convenient. 

3.3. Lattice Simulations 

To simulate the system we discretize space-time and use an adapted leap-frog algo- 
rithm to evolve the system in time (details will be given elsewhere [El)- The 
lattice parameters (all of which are dimensionless) are 

• N_i_, Nj-f the number of lattice sites in the transverse/longitudinal direction 

• g^^a±, a^j, the lattice spacing in the transverse/longitudinal direction 

• Tinit I a± , the time at which the 3-dimcnsional simulations are started 

• (5t, the time stepping size 

• A, the initial size of the rapidity fluctuations 

Of these, only the combinations g^fia^Nx^ = g^^J'L and ajjNjj = (which 
correspond to the simulated system dimensions) have physical meaning; the con- 
tinuum limit is approached by keeping these fixed while sending 6t 0, g^^a± —^ 
0, a,, — > 0. For the 3-dimensional simulations, we still have to choose a value for 
Tinit , which should be such that for A = we stay very close to the result from the 
2-dimensional simulations (for all of which Tinit = 0). Thus, we set 



<F{7j)F{r]') >= A^6{Tj-rj') 



(10) 



Tinit = 0.05 aj.. 



(11) 



but have checked that our results stay the same when choosing Tinit /ci± — 0.025, 0.1, 
respectively. 
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4. Evolution of rapidity-fluctuations and a Weibel instability 
in expanding matter 

An interesting quantity for the longitudinal dynamics is the energy momentum 
tensor component T^^ (see Eq.(|Hll), of which we study the Fourier transform with 
respect to 77, 



where <>_l denotes averaging over the transverse coordinates {x,y). Apart from 
kjj = 0, this quantity would be strictly zero in the boost-invariant (A — 0) case. We 
have checked that this is indeed the case. For finite A, however, T'''' is in general 
non- vanishing for arbitrary fc^ and possesses a maximum amplitude for one specific 
kri- Determining this maximum amplitude for each time step and averaging over 
many random initial conditions, we obtain the curve shown in FigQ The curve is for 
g^fiL = 22.6. In contrast to the curve shown in Ref. [ 6 (for g^fiL — 67.9), where the 
initial seed (violating boost invariance) was chosen to be very small (A ~ 10~^^), 
the seed here is 4 orders of magnitude larger (A ~ 10~^). Nevertheless, for the 
same g^i^L (see Table 1 in Ref. [ 6 -also included below), the fit to the growth rate 
r is consistent with the result quoted in Ref. 

Another feature of our simulations that was not clear from our simulations with 
very small seeds was the flattening of the amplitude, the onset of which is seen in the 
Figure at g^fJ^r w 1000. We have checked (by varying lattice spacing by a factor of 
8) that this phenomenon is rather insensitive to the ultraviolet modes, and appears 
to be a consequence of the "non-Abelianization" of the amplitude. In other words, 
the instability is cut-off when the non-Abelian self interactions of the soft modes 
become important (also seen in Hard-Loop simulations without expansion [ I12j). 
More quantitative studies of this phenomenon will be presented in [113). 

From FigQ] one can see that from g^fir w 150 onwards, there is rapid growth for 
which a best fit (up to times g^/ir 1000) to the functional form cq + ci exp(rfitT'^^) 
gives Ffit = 0.502 ± 0.01 for C3 = 0.5; the coefficients cq, ci are small numbers 
proportional to the initial seed. 

In the presence of a Weibel instability, one expects t^T^^ to grow as exp(~ 7T). 
For a system without expansion, 7 ~ TOoo does not change as a function of r and 
thus the instability manifests itself through modes growing as ~ exp(T). However, 
as argued in Ref. [ (31 , the soft scale moo behaves as ~ l/r in an expanding 
system. Therefore, the functional form of the growth is changed to exp(y^). Our 
results confirm that this functional form is favored by a best fit to our data in Fig^ 

We can confirm this interpretation by also determining moo directly in our 
simulation. This is done by calculating^ the mass gap aj(kj_ = 0) of the gluon 

^The effect of small longitudinal fluctuations on transverse quantities should be rather small. 
Thus, we calculate this mass gap from a 2-|-l-dimensional simulation rather than in the 3+1- 
dimensional case out of computational convenience. 




(12) 
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Fig. 1. The maximum Fomier mode amplitudes of r^T'''' for g^/ii = 22.6, N± = 
Nrj = 32, Lrj = 1.6. Also shown is a best fit with a exp-^r behavior. The flattening 
out of the data at late times is due to the non-Abelian interactions stopping the 
instability growth, 
dispersion relation defined as 



1 / Tr [E.,{k±E,{-k±) + T^E^{k^)E^{-k^)] 
T V Tr [Mk^)A{-k^) + T-2A^(ki)A„(-kx)] ' 



(13) 



which should be proportional to the soft scale a;(k^ = 0) ~ moo- We find [E] 
w(kj_ = 0) = Koy^ g^fi/r, consistent with the expectation from Interpreting 



uj(k± = 0) as the plasmon mass Upi we make use of the relation 



2/3ml m 



to obtain the quantitative estimate moo — '*o-\/3g2/i/(2r) for our simulation. 

If we take the growth rate in the static case and make the change ^statT i{t)t 
with 7(r) = moa{T)/^/2 for the expanding system, we can define the "theoretical" 
growth rate FthcoryVyM^ — '^It. Obtaining Ffit by a best fit to the data (e.g. in 
FigQ for different values of g^^iL, we can compare this result to Fthcory, finding 





Fthcory — Kq 


Ffit 


22.6 


0.526 ±0.003 


0.502 ±0.01 


67.9 


0.447 ±0.003 


0.427 ±0.01 


90.5 


0.49 ± 0.004 


0.46 ±0.04 



However, a consistent treatment would require that the growth rate in the 
expanding case is exp{2 dT'j{T')); if we assume as previously that 7(t') = 
moo(''')/\/2, one obtains an additional factor of 2 in the ratio of T/kq relative to 
Tfit. Understanding this difference of a factor of 2 requires a more careful study of 
the correspondence between the dynamics of expanding fields in our case and that 
in the static HTL case. In particular, it is important to investigate how u>pi that 
we extract from the lattice relates to the plasmon frequency in Hard Thermal Loop 
simulations. Studies in this direction are in progress. 
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